retaildata <- readxl::read_excel("retail.xlsx", skip = 1)
mytimeseries <- ts(retaildata[["A3349873A"]],
frequency=12, start=c(1982,4))
fit <- ets(mytimeseries)
summary(fit)
## ETS(M,A,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.5067
## beta = 1e-04
## gamma = 0.3222
##
## Initial states:
## l = 63.0202
## b = 0.7934
## s=0.9391 0.912 0.9294 1.5281 1.0577 0.9868
## 0.9604 0.941 0.9431 0.901 0.9661 0.9353
##
## sigma: 0.0488
##
## AIC AICc BIC
## 4018.808 4020.494 4085.835
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1576127 13.5932 8.70378 -0.09176148 3.711149 0.4596804
## ACF1
## Training set 0.008030422
autoplot(forecast(fit))
x1 <- window(mytimeseries, end=c(2010,12))
x2 <- window(mytimeseries, start=2011)
f1 <- snaive(x1, h=length(x2))
f2 <- hw(x1, h=length(x2), seasonal='multi')
f3 <- forecast(ets(x1), h=length(x2))
accuracy(f1,x2)
## ME RMSE MAE MPE MAPE MASE
## Training set 7.772973 20.24576 15.95676 4.702754 8.109777 1.000000
## Test set 81.744444 100.00869 82.06667 20.549055 20.669738 5.143067
## ACF1 Theil's U
## Training set 0.7385090 NA
## Test set 0.6830879 1.67023
accuracy(f2,x2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2335309 8.817788 6.393003 0.1048658 3.20614 0.4006455
## Test set 78.8187005 95.779065 78.818701 20.0382885 20.03829 4.9395188
## ACF1 Theil's U
## Training set 0.06611033 NA
## Test set 0.52771581 1.630125
accuracy(f3,x2)
## ME RMSE MAE MPE MAPE
## Training set 0.3404863 8.784116 6.289258 0.1876561 3.128791
## Test set 95.8058083 113.202964 95.805808 24.5580948 24.558095
## MASE ACF1 Theil's U
## Training set 0.3941439 0.01203982 NA
## Test set 6.0040903 0.61158849 1.922081
bicoal %>% ets %>% forecast %>% autoplot
chicken %>% ets %>% forecast %>% autoplot
dole %>% ets %>% forecast %>% autoplot
usdeaths %>% ets %>% forecast %>% autoplot
bricksq %>% ets %>% forecast %>% autoplot
lynx %>% ets %>% forecast %>% autoplot
ibmclose %>% ets %>% forecast %>% autoplot
eggs %>% ets %>% forecast %>% autoplot
ausbeer %>% ets %>% forecast %>% autoplot
autoplot(usnetelec)
No transformation required
autoplot(mcopper)
(lambda <- BoxCox.lambda(mcopper))
## [1] 0.1919047
mcopper %>% BoxCox(lambda=lambda) %>% autoplot
autoplot(enplanements)
(lambda <- BoxCox.lambda(enplanements))
## [1] -0.2269461
# I don't like such strong transformations. Will use 0 instead
enplanements %>% BoxCox(lambda=0) %>% autoplot
autoplot(a10)
(lambda <- BoxCox.lambda(a10))
## [1] 0.1313326
a10 %>% BoxCox(lambda=lambda) %>% autoplot
autoplot(cangas)
(lambda <- BoxCox.lambda(mytimeseries))
## [1] 0.1276369
mytimeseries %>% BoxCox(lambda=lambda) %>% autoplot
f4 <- stlf(x1, lambda=lambda, h=length(x2))
accuracy(f4,x2)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.604219 7.214593 5.176787 -0.2755991 2.525501 0.324426
## Test set 79.914332 97.176605 79.974531 20.2410238 20.263570 5.011954
## ACF1 Theil's U
## Training set 0.02490237 NA
## Test set 0.57946539 1.643056
autoplot(f4) + autolayer(x2, series="Test data")
usnetelec %>% autoplot
usnetelec %>% diff %>% autoplot
usgdp %>% autoplot
usgdp %>% diff %>% autoplot
lambda <- BoxCox.lambda(mcopper)
mcopper %>% BoxCox(lambda=lambda) %>% autoplot
mcopper %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot
enplanements %>% log %>% autoplot
enplanements %>% log %>% diff(lag=12) %>% autoplot
enplanements %>% log %>% diff(lag=12) %>% diff %>% autoplot
visitors %>% autoplot
lambda <- BoxCox.lambda(visitors)
visitors %>% BoxCox(lambda=lambda) %>% autoplot
visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot
visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% diff %>% autoplot
mytimeseries %>%
BoxCox(lambda=0.12) %>%
diff(lag=12) %>%
diff(lag=1) %>%
autoplot
wmurders %>% autoplot
wmurders %>% log %>% autoplot
fit <- auto.arima(wmurders, lambda=0, stepwise=FALSE, approximation=FALSE, max.p=10, max.q=10, max.order=10)
checkresiduals(fit)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,2,1)
## Q* = 10.562, df = 8, p-value = 0.2278
##
## Model df: 2. Total lags used: 10
forecast(fit) %>% autoplot
wmurders %>% ets %>% forecast %>% autoplot
wmurders %>% ets(lambda=0, model="AAN") %>% forecast %>% autoplot
lambda <- BoxCox.lambda(mytimeseries)
auto.arima(mytimeseries, lambda=lambda)
## Series: mytimeseries
## ARIMA(1,1,1)(0,0,2)[12] with drift
## Box Cox transformation: lambda= 0.1276369
##
## Coefficients:
## ar1 ma1 sma1 sma2 drift
## 0.5676 -0.9565 0.9543 0.5450 0.0096
## s.e. 0.0523 0.0183 0.0528 0.0357 0.0022
##
## sigma^2 estimated as 0.02727: log likelihood=140.48
## AIC=-268.95 AICc=-268.73 BIC=-245.31
(arimamod <- auto.arima(mytimeseries,
D=1, lambda=lambda,
stepwise=FALSE,
approximation=FALSE))
## Series: mytimeseries
## ARIMA(1,1,2)(0,1,1)[12]
## Box Cox transformation: lambda= 0.1276369
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.4147 0.0703 -0.3325 -0.5811
## s.e. 0.1898 0.1822 0.0725 0.0555
##
## sigma^2 estimated as 0.009902: log likelihood=326.42
## AIC=-642.84 AICc=-642.67 BIC=-623.3
newretaildata <- readxl::read_excel("8501011.xls", sheet=2, skip=9)
mynewdata <- ts(newretaildata[["A3349873A"]],
frequency=12, start=c(1982,4))
mynewdata <- window(mynewdata, start=2014)
lambda <- BoxCox.lambda(mytimeseries)
(etsmod <- ets(mytimeseries))
## ETS(M,A,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.5067
## beta = 1e-04
## gamma = 0.3222
##
## Initial states:
## l = 63.0202
## b = 0.7934
## s=0.9391 0.912 0.9294 1.5281 1.0577 0.9868
## 0.9604 0.941 0.9431 0.901 0.9661 0.9353
##
## sigma: 0.0488
##
## AIC AICc BIC
## 4018.808 4020.494 4085.835
f1 <- snaive(mytimeseries, h=length(mynewdata))
f2 <- hw(mytimeseries, h=length(mynewdata), seasonal='multi')
f3 <- forecast(etsmod, h=length(mynewdata))
f4 <- stlf(mytimeseries, lambda=lambda, h=length(mynewdata))
f5 <- forecast(arimamod, h=length(mynewdata))
checkresiduals(f1)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 804.66, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
checkresiduals(f2)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 40.245, df = 8, p-value = 2.884e-06
##
## Model df: 16. Total lags used: 24
checkresiduals(f3)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,M)
## Q* = 26.597, df = 8, p-value = 0.0008296
##
## Model df: 16. Total lags used: 24
checkresiduals(f4)
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(M,A,N)
## Q* = 69.123, df = 20, p-value = 2.531e-07
##
## Model df: 4. Total lags used: 24
checkresiduals(f5)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)(0,1,1)[12]
## Q* = 16.502, df = 20, p-value = 0.685
##
## Model df: 4. Total lags used: 24
accuracy(f1,mynewdata)["Test set","RMSE"]
## [1] 77.30981
accuracy(f2,mynewdata)["Test set","RMSE"]
## [1] 52.98594
accuracy(f3,mynewdata)["Test set","RMSE"]
## [1] 53.71773
accuracy(f4,mynewdata)["Test set","RMSE"]
## [1] 39.35806
accuracy(f5,mynewdata)["Test set","RMSE"]
## [1] 29.27015
autoplot(f5) +
autolayer(mynewdata, series="New data")